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I. INTRODUCTION 



Understanding complex nonlinear phenomena in magnetized plasmas increasingly relies 
on the use of numerical simulation as an enabling tool. The development of a robust pre- 
dictive capability requires numerical models which are verified through comparison with 
analytic calculation and validated through comparison with experiment 1 . A tractable ana- 
lytic problem useful for verification of numerical models of plasma turbulence and transport 
is linear stability 2 " 4 -. Understanding of linear instabilities in a set of model equations forms a 
framework for developing physical insights and mathematical apparatus that can be further 
used for attacking a more difficult nonlinear problem. 

This paper presents a study of linear gradient- driven instabilities in a cylindrical mag- 
netized plasma using the Braginskii two-fluid model. This work was undertaken with two 
motivations: (1) to gain understanding of the character of linear instabilities in the Large 
Plasma Device (LAPD) at UCLA 5 and (2) to verify linear calculations using the BOUT 3D 
Braginskii fluid turbulence code 6 in cylindrical geometry. The BOUT code was originally 
developed in the late 1990s for modeling tokamak edge plasmas; the version of the code used 
in this study is described in detail by Umansky-. 

Instability drive in LAPD comes from plasma pressure gradients^ and strong azimuthal 
flow which can be externally driven through biasing*^. These free energy sources can 
drive resistive drift waves^ 1 - and Kelvin-Helmholtz and rotational interchange instabilities^. 
The Kelvin-Helmholtz instability and unstable drift-Alfven waves have been experimentally 
observed in LAPD^ 1 ^ 1 ^ 4 -. The analytic calculations and verification runs on BOUT are 
performed using LAPD-like profiles of plasma density, temperature and plasma potential 
(and therefore cross-field E x B flow). It is found that all three modes (drift waves, Kelvin- 
Helmholtz, rotational interchange) can be important in LAPD plasmas. Detailed comparison 
with solutions of the analytic dispersion relation demonstrates that BOUT accurately re- 
produces all characteristics of linear modes in this system. This work forms the foundation 
for nonlinear modeling of turbulence and transport in LAPD, initial results of which will be 
presented in a companion paper—. 

This paper is organized as follows. Section [III introduces the LAPD geometry and presents 
the fluid model used for calculations of linear instabilities. Section II I II discusses the imple- 
mentation of these equations in the BOUT code, including a discussion of techniques used 
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to extract characteristics of linear instabilities. Comparison of BOUT calculations to ana- 
lytic linear eigenmode solutions are presented in Section [IV] for three instabilities: resistive 
drift waves, Kelvin- Helmholtz and rotational interchange modes. Section |V] discusses the 
linear stability of experimentally measured LAPD profiles against these three instabilities 
and a discussion of the similarity to experimental observations. The effect of ion-neutral 
collisions on the linear solutions is discussed in Section IVII A summary of the paper is 
presented in Section IVHI Appendices are provided which cover: a derivation of the specific 
set of fluid equations used in this work (Appendix |A|) ; a derivation of the vorticity equation 
used in BOUT (Appendix [B]) ; and a list of parameters and boundary conditions used in the 
verification study (Appendix 0. 

II. GEOMETRY AND PHYSICS MODEL 

The geometry used in this study is that of the LAPD: a ~ 17 m long cylindrical magne- 
tized plasma with typical plasma radius (half- width at half-maximum) of a ~ 30 cm (vacuum 
chamber radius r = 50 cm). Typical plasma parameters in LAPD for a 1 kG magnetic field 
are shown in Table [H 

The configuration is modeled as a cylindrical annulus to avoid the singularity of cylindrical 
coordinates near the axis in the BOUT numerical implementation (Fig [I]). Using the scheme 
shown in Fig [U LAPD geometry can be completely described within BOUT framework 
without major modification of the core code. The only change related to geometry in the 
code that is necessary is the implementation of the full cylindrical Laplacian operator to 
extend the simulation domain closer to the magnetic axis. 

The magnetic field is taken uniform, directed along the cylinder axis. The axial boundary 
conditions are taken periodic for simplicity. A more realistic model should include the end- 
plate sheath boundary conditions, supporting potentially important wall-driven instabilities; 
this will be the subject of future work. Radial boundary conditions used here are either zero 
value or zero radial gradient. 

A Braginskii two-fluid model 16 is used in the analytic and BOUT calculations for insta- 
bilities in LAPD. As evident from Table [IJ collisions are important in LAPD plasmas: the 
electron collision mean free path is much smaller than the system size parallel to the mag- 
netic field, A ei <C L\\. Therefore for long parallel wavelength, low frequency modes (a; <C Qi) 
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FIG. 1. (Color Online) Schematic view of LAPD geometry representation in the BOUT code. 
The poloidal direction of the tokamak geometry becomes the axial direction z in LAPD, and the 
toroidal coordinate of a tokamak corresponds to the azimuthal angle 9 in LAPD. 



Species 


4 He 


fci 


380 kHz 


Z 


1 


Pi 


0.2 cm 


n 


2.5 x 10 12 cm~ 3 


Ps 


0.5 cm 


T e 


5 eV 




7.4 x 10 6 1/s 


Ti 


<leV 


fa 


5 x 10 5 1/s 


Bo 


0.1 T 




1.2 x 10 3 1/s 


L \\ 


17 m 




13 cm 


a 


~ 0.3 m 




~ 4 x 10 4 rad/s 



TABLE I. Typical LAPD parameters 



considered here, it could be argued that the use of a collisional fluid theory is justified. How- 
ever, it should be noted that the quantity most important for evaluating the importance 
of kinetic effects is the ratio of the parallel wave phase speed to the thermal speed of the 
particles, and for drift-type modes and Alfven waves in LAPD this can be near unity for the 
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electrons^ 1 ^. Strong collisions can disrupt velocity-space resonant processes and it might 
be expected that a fluid description becomes accurate even for v$ ~ v t h as k\\X e i —> 0, as 
has been shown for ion acoustic waves through Fokker-Planck calculations^. The present 
work is part of an ongoing effort to evaluate the validity of a fluid model (in particular 
that implemented in BOUT) in describing turbulence in LAPD. A goal of this study is to 
determine whether (and how) fluid simulations can fail to describe plasma behavior, and 
kinetic effects are likely to delineate when failure occurs. 

The fluid equations used here represent conservation of density, electron and ion momen- 
tum and charge: 



where p e = nk^T e . A friction term due to ion- neutral collisions (elastic and charge-exchange) 
is included in the ion momentum equation. All terms involving finite ion temperature effects 
are neglected. The friction forces in the electron momentum equation are due to electron-ion 
(v e i) and electron- neutral collisions (v en ). However, as Coulomb collisions are dominant for 
the electrons [y ei ^> u en ), electron-neutral collisions are ignored. 

The following simplifying assumptions are made, which are relevant for LAPD plasma 
parameters: constant magnetic field B = B z, v\\ e ^> T e ^> T iy and no background par- 
allel flows. In addition, it is assumed that the instabilities do not generate perturbations in 
the electron temperature. Throughout the paper plasma density, temperature and magnetic 
field are normalized to reference values n x , T ex (chosen as the maximum of the correspond- 
ing equilibrium profiles), and Bo, the axial magnetic field. Frequencies and time derivatives 
are normalized to Q{ X = eBo/rriiC: dt = dt/fl{ X , w = u/Qi X ; velocities are normalized to 
the ion sound speed C sx = sjT^fm^ lengths - to the ion sound gyroradius p sx = C sx /Vt ix ] 
electrostatic potential to the reference electron temperature: = e<p/T ex . Further the " " " 
symbol for dimensionless quantities will be dropped for brevity of notation. 



(dt + v e • V) n 







(1) 




(2) 
(3) 
(4) 
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Combining Eqs.([T]|4]) and linearizing (see Appendix |A]), one obtains: 

d t N + b x V±0 O • ViV = -b x V±0 • VAo - iV V[|Uj[ a 

d t v\\ e + b x V±<f) ■ Vv\\ e = -/i-^-V||iV + /iV||0 - u e v\\ e (5) 

iv 

AW Flic = -V ± • (N Q d t V±(f> + d t NV x <Po 

+b x V ± O • V (AW±0o) 
+b x V ± O • V (N V ± <f>) 
+b x V ± • V (N V ± <f> ) 
+b x V ± 0o • V (iW ± o ) 

+A 7 " O Z/i„V_L0o + N V in V±(f) + Nu in V±(p J , 

where Ao, 4>o, T e0 are zero-order (equilibrium) quantities and N, <fi, v\\ e , are first order pertur- 
bations; n = mi/m e . 

Note that Eqs. (J5j) contain a zero-order term, Vj_ ■ (^mAoVj^o), which restricts the choice 
of background profiles in the presence of neutrals. If this term is not zero for a particular 
choice of Nq(t) and </>o( r ) functions, then the plasma is not in mechanical equilibrium. In 
such a case, an extra zero-order force, e.g., from externally applied radial electric field, should 
be added to the momentum equation to balance the force of friction with the neutrals that 
slows down the plasma rotation. 

Next, Eqs. (jSJ) are projected onto cylindrical coordinate system (r,9,z). Solutions of 
the form /(x) = f (r) exp(im,g9 + ik\\z — iuot) are sought, where fcn = 2im z /L\\, n z is the 
parallel mode number. Denoting /' = d r f and introducing the Doppler-shifted frequency 
Co = u — ^r4>' , the ID equation for radial eigenfunctions of the perturbed potential <f>(r) can 
be written: 

C 2 (r)0" + C 1 (r)0 / + C o (r)0 = O, (6) 

where the coefficients Cj(r) are functions of equilibrium quantities and of Co (full expressions 
for Ci are presented in Appendix |A|) . 

Equation (jBJ) is a 2nd order ordinary differential equation (ODE) in r. Supplemented 
with proper boundary conditions on the radial boundaries it forms a well-posed eigenvalue 
problem. In general Eq. ([6]) has to be solved numerically, due to the complex form of the 
coefficients Cj. Note that although Eq. ([6]) is better suited for theoretical analysis than the 
original system, Eqs. (jSJ), a complication for practical numerical solution of Eq. ([6]) is that 

6 



the eigenvalue u enters nonlinearly the coefficients Cj. Therefore, a numerical solution for 
the eigenvalues is easier to carry out using the original system Eqs. (jSJ), which can be cast to 
a standard linear algebra eigenvalue problem amenable to solution by a standard eigenvalue 
package. 

III. SOLVING BY TIME-EVOLUTION WITH BOUT 

The present version of the BOUT code^ is a rather general framework suitable for inte- 
gration of a system of time-evolution PDEs in 3D space of the form d t f = F(f, x), where 
the right-hand-side F contains a combination of spatial differential operators applied to the 
state vector f . The right- hand- side F is discretized on a spatial mesh by finite-differencing, 
which results in a system of ODEs that are integrated in time by an implicit ODE solver 
package PVODE^. 

For the calculations presented here, the following set of equations are used in BOUT 
which are equivalent to Eqs. (jTEj): 

d t N = -v E - VN-V\\(v h N) (7) 

T 

d t V\\ e = -Vjs ■ VV|| e - fl-^V\\N + A*V||0 - V e V\\ e (8) 

iVo 

d t w = -v E ■ Vro - V||(iVu|| e ) + b x VN • Vu|/2 - v in w (9) 
where the potential vorticity 

w d = f v ± • (A^V ± 0) (10) 

is introduced. While the variables N, v» e and w are advanced in time, Eq. ( ITU]) is inverted 
on each evaluation of the right-hand side of Eqs. (JTJIH]) to reconstruct the perturbed potential 
(f) from vj. 

The vorticity evolution equation, Eq. (JUJ), replaces the current continuity equation (j4]) 
in BOUT. Derivation of this form of the vorticity equation from Eq. (j4]) is presented in 
Appendix [B] Note that Eq. ([UJ is equivalent to expression (76) in Simakov and Catto^S, 
apart from the ion-neutral collision term which is included in this work and all the terms 
involving ion temperature which are neglected in the present work. The third term in 
the right-hand side of Eq. (]9j is important even in the linear regime if strong background 
flows are present. Thus, it is essential in both linear and nonlinear simulations of LAPD 
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experiments, which typically have strong azimuthal flows on the order of Mach number M 
~ 0.2 for spontaneous flows and M ~ 1 in bias-induced rotation experiments^. 

To compare BOUT solution of the initial-value problem with the direct solution of eigen- 
value problem corresponding to the discretized dispersion relation (jSJ), the equations are 
linearized (Eqs. [7 HT0]) and advanced in time using BOUT from a small initial seed pertur- 
bation. BOUT produces perturbations, in this case, of density and vorticity/potential, as 
functions of 3D space and time. 

A specific azimuthal mode number nig is selected by Fourier filtering in the azimuthal 
angle during the BOUT simulation. The parallel wave number k\\ is set by the length of the 
device and the periodic boundary conditions in the parallel direction. 

The radial form of the numerical solution is dominated by the fastest growing radial 
eigenmode. Once the solution "locks in" to the fastest mode, we calculate the growth rate 
by fitting the time evolution of the volume-averaged amplitude of potential fluctuations to an 
exponential. The frequency of the mode is then calculated by fitting the perturbed potential 
(with the exponential growth factored out) with a sine wave at each spatial position. 

IV. VERIFICATION OF BOUT AGAINST EIGENVALUE SOLUTION 
A. Electrostatic resistive drift wave 

In the absence of strong flows, the resistive drift mode is likely to be the primary instability 
in LAPD. In this section, the BOUT solution in LAPD geometry is verified using a reduced 
subset of fluid equations Eqs. (17HT01 which support only the resistive drift instability branch. 

The simplest model of the resistive drift wave can be written as a subset of the system 
Eqs. (EU): 

d t N = -v E ■ ViVo 
T o 

d t V\\ e = — jU— ^V||A^ + yUV||<^ — U e V\\ e (11) 
d t VJ = -V||(JVV|| e ) 

These equations can be combined together to form a well-known local dispersion relation^- 
that assumes ID dependence of the background density with constant gradient length L n = 
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N (x)/N^x): 




w S\ . a \\ ( w 
1 »-+ - 



) 



2 



(12) 



where 



BOUT calculations were first verified on this simple local solution, Eq. ( TT2l) . finding good 
agreement for a range of plasma parameters. Due to its simplicity, this solution provides 
useful insight into the behavior of the growth rates and frequencies. In a bounded plasma, 
the dispersion relation Eq. (|12p together with a set of boundary conditions yields a set of 
discrete linearly unstable modes. Among these discrete modes, the fastest growing one is 
the mode that corresponds to dimensionless parameter cry /a;* closest to 1. 

The next, and more interesting, step is to compare BOUT calculations to the eigenvalue 
solution of the full non-local drift wave problem Eqs. ([Ull]). Here all terms in Eqs. fl71 fTU|) 
are retained. There is no background potential </>o for these calculations which eliminates 
the Kelvin-Helmholtz and the rotation-driven interchange instabilities, and only allows for 
the drift wave solution. There is no simple analytic dispersion relation in this case. For 
comparison, we are using the direct numerical solution of the linear problem Eq. (JH]) obtained 
with an eigenvalue solver, as described in Section HI] The results of this comparisons for 
cylindrical geometry with relevant to LAPD parameters and profiles are presented in Fig. (J2J). 
BOUT recovers the frequencies and growth rates for a range of magnetic field values (B = 
0.04, 0.08, and 0.12T). There is one-to-one correspondence between the eigenvalues found 
by the analytic solver and the BOUT solution. Typically, the discrepancy between the two 
methods is less than 2% for radial grids of 50 points, and the results converge with grid size. 
For comparison, frequencies and growth rates for longer wavelengths, n z = 0.5 (fundamental 
mode), are also shown (dashed lines). 

As an example of numerical convergence, the relative error of the growth rate and fre- 
quency as a function of radial grid size h = 1/N r is shown in Fig. |3j The relative error is 
defined here as the difference between BOUT solution and the projected value at h = 0, 
<5 7 = 1 7 — r yh=o\/lh=o, and analogously for the frequency uj. The growth rate and frequency 
extracted from the initial-value simulation converge approximately quadratically in h. The 
difference between the BOUT solution at h — > and the eigenvalue solver result is 0.43% 
for the frequency and 0.23% for the growth rate. This residual error is due to the lim- 
ited numerical resolution in the azimuthal and parallel directions (both remain fixed at 16 
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FIG. 2. (Color Online) Frequencies (left) and growth rates (right) of the resistive drift wave in 
LAPD configuration for experimental density profile, without equilibrium flows. Analytic results 
are shown with lines for B = 0.04 T (solid), 0.08 T (dash), 0.12 T (dotted); corresponding BOUT 
results are shown with squares, circles and triangles. The two groups of lines correspond to the 
axial harmonics n z = 1 and n z = 0.5. mg values on top axis are given for B = 0.04T. 

grid points for this convergence study) and slight differences in the representation of the 
underlying equilibrium in BOUT and the eigenvalue solver. 

B. Kelvin-Helmholtz instability 

LAPD plasmas often involve large azimuthal flows, especially with biasing of the vacuum 
vessel wau^ 1 ^. The flows in the experiments with externally applied radial bias can reach 
Mach number of about 1, or vg ~ 10 6 cm/s. These speeds are much higher than the 
typical phase velocity of the drift wave, Vd ~ 0.5 x 10 4 cm/s. Also, the growth rates of the 
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FIG. 3. (Color Online) Relative error of the growth rate and frequency determined from the 
initial-value simulation as a function of the radial grid size, indicating 2nd order convergence. 

instabilities generated by bias-induced flows can be comparable to that of the drift wave 
(see Section IVj, therefore it is essential to include these flows in the model. 

Kelvin-Helmholtz (KH) instability, driven by sheared flows, represents an interesting 
case for BOUT simulations in LAPD geometry, and provides a test of the implementation 
of the terms involving 0o in BOUT. Observations of KH instability in LAPD plasmas have 
previously been reported by Horton et al.— . 

A simple model for the KH instability can be obtained from the charge conservation 
equation, Eq. fljj). Assuming no variation of equilibrium or perturbed quantities along the 
magnetic field (flute modes), only the polarization current contributes to this equation: 

J± = e,(v a - v«J = = (9, + • V) V ±4 > (13) 

For simplicity, the case of constant plasma density and magnetic field is considered. The 
charge conservation equation can then be written as 

(«9i + VjE - V)Vl0 = O (14) 
11 



Linearizing Eq. (|T4|) for slab geometry with periodic coordinate y we obtain the eigenvalue 
equation*^: 

where the solution is assumed of the form 0(r, t) = 0(x) exp(ik y y — iut) 

Analytic solution of this equation can be found for a specific choice of stream function 
0o by matching and its derivative jump at the points of singularity: 

' 0, x<-l 

I x 2 /2 + x+l/2, -l<x<0 
<> (x) = { (16) 
-x 2 /2 + x + 1/2, 0<a;<l 

1, x > 1 

For direct comparison with BOUT, a solution must be found with boundary conditions 
imposed on a finite interval. We consider boundary conditions 0(— 2) = 0(2) = const. In 
this case, the eigenvalues are 

u = (e 2k y - l)/(2 + 2e 2k y) (17) 
for the neutrally stable branch, and 



(l - e Ak * + 2k y + 2 k y e Aky ± ^G(k y )j /(4 + 4e 4fca 



G(k y ) = 9 - 16e^ + 14e 4fc » - 16e b/Cy + 9e 8/Cy + 12^ - I2k y e* k * 
+4k 2 y + 8k 2 y e 4k » +Ak 2 y e 8k y. 

for the stable/unstable branches. One of the branches is unstable for < k y < 1.815, 
maximum growth rate is 0.2346 at k y 1.241. This result is similar to the calculation 
presented by Horton et al. for slightly different boundary conditions 9 . 

This instability is found with BOUT by solving Eq. (fl4|) written in terms of vorticity: 

d t vo + v E -Ww = (18) 

Eq. (fT8"|) is explicitly linearized in BOUT and solved in slab geometry with the same bound- 
ary conditions </>(— 2) = 0(2) = const. In BOUT, slab geometry is approximated as a small 
azimuthal segment of a large aspect ratio thin annulus. The exponential growth rate and 
the mode frequency is extracted from the time evolution of the perturbed potential 0. Using 
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FIG. 4. (Color Online) Analytic solution (frequency - dashed line, growth rate - solid line) and 
BOUT simulations (circles) for the Kelvin-Helmholtz instability in slab geometry, (fro profile for 
this case is given by Eq. (fT6|) . 

this method, the frequencies and growth rates of the direct eigenvalue solution are recovered, 
as shown in Fig. HI 

Note that the third derivative of <fio that is present in Eq. (jl~8l) (as can be seen from 
Eq. ( 1T5|) ) is singular, but it does not directly enter BOUT equations. BOUT uses 4>o profile 
as input, which is a smooth function. Therefore, the code has no difficulty reproducing the 
analytic solution even though it implies a singularity in the 0q" profile. 

Next, to make a calculation relevant to the experiment, the KH instability in LAPD 
geometry is considered using the experimental density profile and a model 0o( r ) profile with 
amplitude values relevant to the experiment. The background potential profile is similar to 
expression (fl~6l) . but the delta-functions in are replaced by Gaussians (exact expression 
is given in the Appendix EH Eq. ( IClj) ). This calculation represents a strong test of the terms 
involving background flows in Eqs. (l9l fT0l) since some of these terms only contribute when 
both VN and V0o exist. Note also that with non-constant N (r), the density perturbation 
is not zero, unlike in the situation considered above. There is no analytic solution in this 
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case, therefore we compare the BOUT solution with the results of the eigenvalue solver for 
the system of equations (fUHJ). The comparison is presented in Fig. [5j The result is similar 
to the previous KH case, with a cutoff in perpendicular wavenumber. In LAPD geometry, 
for this particular choice of profiles, this cutoff translates into itlq w 8; the KH mode is 
stable above this value. BOUT reproduces the direct eigenvalue solution with a very good 
accuracy of < 2% for a 100 point radial grid size. 




O06 0\08 010 012 0J4 016 0.18 0.20 
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FIG. 5. (Color Online) Eigenvalues (frequency - dashed line, growth rate - solid line) of the 
Kelvin-Helmholtz instability as a function of perpendicular wavenumber. Circles - BOUT results. 
Cylindrical geometry, experimental density profile, LAPD plasma parameters. 

C. Interchange instability 

Strong azimuthal flows in LAPD not only affect the frequency of the waves through a 
Doppler shift, but can also modify the growth rate even for uniform rotation due to the 
induced centrifugal force. In this section the rot at ion- driven interchange mode is considered 
in the presence of background density gradient. 

To separate the interchange mode from the other instabilities, the parallel wave number 
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is set to zero (this removes the drift wave branch of the dispersion relation) and a uniform 
rotation profile is chosen with normalized rotation frequency Q: <po = (px^/^max = ^ r2 /2 
(this removes the KH mode). Assuming an exponential equilibrium density profile with 
gradient scale length L n , Eq. fl^]) is written as 



m 



2 



+ - —<// - — + — U - — <f> = 0, (19) 
r L n rL n \oo z lo J r A 

where d is the Doppler shift, d = mcf)' Q /r = 2mcp x /r' 2 nax = mfl, lo = uj — d, and L n = N /Nq. 

There are two tractable limits where analytic solution can be found, kL n 3> 1 (slow 
variation), and kL n <^ 1 (sharp interface). 

For small density gradient, the 1/L n term can be dropped when compared with the 1/r 
term. Employing a change of variable x = y/r, Eq. f lT9|) is rewritten as Bessel's equation: 

x 2 <\)\x) + x4>'{x) + {AC 2 - Am 2 x 2 ) <f>[x) = (20) 
1 f d 2 2d\ 

where C 2 = ——— I + — ] . The solution is given as a sum of Bessel functions of the first 



L n \oo 2 lo J 
and the second kind: 

(f>(x) = d J 2m (-2Cx) + C 2 Y 2m (-2Cx). (21) 

The dispersion relation is obtained by imposing the boundary conditions 0(r m j n ) = (p(r max ) = 
on this function. For simplicity r min = is chosen. Y 2m (x) diverges at the axis, so the 
dispersion relation in this case is given by the condition 

J 2m (-2C V ^~) = 0. (22) 

For large mg the position of the first zero of the Bessel function J m {x) can be estimated 24 
as mg (e.g. j m = 36.1 for m = 30 and the relative error monotonically decreases for larger 
mg). This results in a simple approximate equation for the interchange eigenmode: 

— Cy/ r max = m, (23) 

which yields the approximate dispersion relation (again using L n 3> r max ) 



u = mQ±tQj r -^ (24) 

Note that the growth rate 7 = Vtyjr max / L n can be obtained from the well known disper- 
sion relation of the Rayleigh- Taylor instability driven by gravity, 7 ~ a/ 'g/L n , if gravity is 
replaced by the centrifugal force of the rotation, g = Q 2 r max . 
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FIG. 6. (Color Online) Interchange mode in a rotating cylinder for the case of exponential density 
profile (a) and piecewise-constant density (b). Solid line - eigenvalue solution, dashed line - 
asymptotic at large k±, dashdot line - exact analytic solution for case (b), circles - BOUT results. 

The growth rate given by Eq. (I24p is independent of and represents an asymptotic 
solution for large m$. This asymptotic solution and the exact solution of Eq. (fT9l) are shown 
in Fig. E](a). 

Another limit where a simple analytic solution can be found is the case of a piecewise- 
constant density profile with a sharp interface, A (r) = Aq for r < r and N (r) = A 2 for 
r > r , r = r max /2, Aq > N 2 . Eq. (TT91) at r ^ r then becomes 

<P" + V - ^ = 0, (25) 

with the general solution 0(r) ~ r ±m . Matching the values of 0(r) at the interface, applying 
the boundary condition at the conducting shell 4>{r max ) = 0, and integrating Eq. (fT9|) in a 
small region near the interface to account for the jump in A and 0', we obtain the dispersion 
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relation 

Co = Q (^—a ± V a 2 — mc^j , (26) 

where a = A(2 2m - l)/(2 2m - A) and A is the Atwood number A = (N 1 - N 2 )/(N 1 + N 2 ). 

In the limit of large m# cylindrical effects become insignificant and the growth rate con- 
verges to that of the gravity-driven Rayleigh- Taylor instability in a slab for two fluids with 
sharp interface, 7 ~ \/VL 2 mA = ^kgA. The exact solution and the asymptotic solution at 
large irtg are shown in Fig. E^b). The solid line represents the eigenvalue solution of the 
system ([5]) where the piecewise-constant density profile is approximated with tanh function. 
At higher mg numbers (here, at m > 20, or kgp s > 1) the finite width of the interface 
region becomes important compared to 1/ke, so the numerical (eigenvalue) solution starts 
to deviate from the analytic solution (|2"B1) . 




FIG. 7. (Color Online) Interchange mode (/cy = 0) destabilized by uniform rotation and drift- 
interchange mode (/c|| = 2n/L). 

The system of time-evolution equations used in BOUT to reproduce the interchange mode 
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can be obtained from Eqs. (JTHnj) by setting k\\ = 0: 



d t N 



d t w 



v E ■ VN 



v E ■ Vro + b x VA^ • Vv E /2 - v in w 



(27) 
(28) 



where all variables (N, v#, vo) contain both the equilibrium part and the fluctuating com- 
ponent. These equations are linearized in BOUT, and solved with the same parameters as 
used in the two analytic examples discussed above (Fig. |6]). BOUT simulation recovers the 
interchange mode solution for both limits (slowly varying exponential and piecewise constant 
profiles) which verifies the correct implementation of the new terms involving background 
flows in BOUT. 

In order to investigate the effect of uniform rotation on the interchange and drift- 
interchange instabilities in LAPD plasmas, a configuration with experimental density profile 
and 4>o(r) ~ r 2 is considered. The results of this calculation are presented in Fig. (jJJ) as 
a function of rotation velocity. Two axial harmonics are shown: n z = (pure interchange 
mode) and n z — 1 (drift-interchange instability). At <po = 0V, n z — 1 branch corresponds 
to a pure drift mode. As O increases, the frequency of this mode is Doppler-shifted and the 
growth rate is modified by the centrifugal force. At large rotation velocities, the disparity 
between the large real part of the frequency and the small growth rate is hard to resolve 
numerically using an initial- value code, so the BOUT results slightly deviate from the direct 
eigenvalue solution of the dispersion relation. 

V. LINEAR INSTABILITIES IN LAPD 

Now that simple analytic solutions for each of the instabilities supported by Eqs. flTJHJ) 
have been presented, linear instabilities for LAPD parameters and experimental profiles will 
be considered and the growth rates for the different mode branches will be compared. 

Fig. (jEJ) shows the growth rates and frequencies of the KH, drift and interchange modes 
for LAPD parameters using experimentally measured density profiles. The complete set of 
parameter values including the polynomial fit of the experimental density profile is shown 
in Appendix O Three different model background potential profiles are chosen here to 
separate the instability branches: same profile as used in Fig. [5] for the KH mode (given by 
Eq. ( 1C1|) ). uniform rotation profile <pa{r) ~ r 2 for the interchange mode, and zero potential 
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FIG. 8. (Color Online) Kelvin-Helmholtz, drift and interchange branches of the dispersion relation 
for LAPD parameters as a function of azimuthal mode number. Left: equilibrium profiles, right: 
frequency and growth rate of the instability. For KH and IC n z =0, for DW n z =0.5. For DW case 
00 ( r ) = 0, for IC instability (f>o(r) ~ r 2 (uniform rotation), for KH mode 4>o(r) is given by Eq. (|Clj) . 
As a reference, the experimentally measured 4>lapd profile is shown in dotted line (left). 

for the drift wave instability. The magnitude of the radial potential drop in the KH and IC 
cases is of the same order as the measured value in biased discharge experiments^. Even 
though a direct comparison of the three solutions is not possible because the background 
flow profiles and axial mode numbers are not the same, it is still informative to note that 
the growth rates of all three branches of instability are of similar magnitude. Therefore, all 
three instabilities can potentially compete in LAPD plasmas. 

Similar results are observed in a calculation of linear growth rates when using self- 
consistent, experimentally measured profiles of density, electron temperature and flow 
(Fig. [9j). Two cases are considered here, biased and unbiased plasma discharges (with 
and without bias-driven azimuthal edge flow}^. In the unbiased configuration, the az- 
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FIG. 9. (Color Online) Equilibrium profiles of the density, electron temperature, potential and 
Mach number (left) and the fastest growth rates (right) of the perturbation with axial mode 
number n z = 0,0.5,1. Top: unbiased LAPD discharge. Bottom: LAPD discharge with applied 
radial bias. Experimental data taken from Maggs et al.— . 



imuthal flow values are much smaller than in the biased case, so we use zero azimuthal flow 
for this calculation. In the unbiased case (Fig. [91 top), only the drift wave branch is present, 
with comparable maximum growth rates for n z — 1 and n z = 0.5. In the biased case (Fig. EjJ 
bottom), the growth rates at nig < 10 for the three harmonics n z = 0, n z = 0.5 and n z = 1 
are comparable. From the eigenfunction analysis, it can be concluded that n z = 1 harmonic 
is predominantly interchange at mg < 5, then drift wave- like at 5 < rriQ < 17 and again 
IC-like at higher m$. An example of the eigenfunctions of the potential perturbation for the 
biased case is shown in Fig. [TUJ At itlq = 3 and nig = 20, the axial mode n z = 1 is localized 
near the edge of the plasma where the azimuthal flows are strongest (see Mg profile in 
Fig. El bottom), which is consistent with the rotational interchange instability. At nig = 12, 
the n z = 1 harmonic eigenfunction is localized near r ~ 28 cm, where the gradients of the 
density and electron temperature are strongest, which indicates the drift- wave- like character 
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of the mode. 

The real frequencies of these modes are consistent with experimental observation; in the 
unbiased case, at the peak of the growth rate for n z = 0.5 (m ~ 20) the mode frequency 
is / = 4.7kHz which is in the heart of the measured broadband fluctuation spectrum in 
unbiased plasmas, although a lower mode number would be consistent with the measured 
correlation function^. In the biased case, the local maximum of the growth rate of n z = 0.5 
mode is at me = 8, which is also consistent with measured LAPD value of me < 10 10 . 
Higher growth rates at large me might not be relevant when viscosity effects are included 
in the calculation, since high k± modes will be damped by viscosity. The computed linear 
eigenf unctions are consistent with the observed fluctuation profiles in the unbiased case, 
localized to the density gradient region. In the biased case, eigenf unctions localized to the 
region of strong density gradient are found as well as flow-driven modes that are localized 
to the far edge away from the strong gradient region. The latter is consistent with the 
observation of increased electric field fluctuations in the far-edge plasma with increased bias 
(see Fig. 9c of Ref. 10). The linear prediction that Kelvin- Helmholtz and/or rotational 
interchange might be the dominant instabilities in the biased case could also be consistent 
with measurements of the cross phase between density and electric field fluctuations. In 
going from unbiased to biased plasmas in LAPD, a dramatic change in the cross phase is 
observed, which could be consistent with a change in the dominant instability 10 . 

A detailed comparison with the experimental data requires nonlinear analysis and simu- 
lation, which is the subject of a companion paper—. However, it is still illustrative to apply 
quasilinear theory or mixing length arguments^ 1 ^ using the linear calculation results pre- 
sented above. For drift waves driven by the background density gradient the mixing length 
estimate assumes that the saturation is reached when the perturbed gradients become com- 
parable to the equilibrium gradients: 



nk± ~ n /L, 



'n 



(29) 



so that 



n/n ~ q4>/T e ~ 1 



fk±L. 



(30) 



The weak turbulence theory modifies this estimate by a factor 




n 



/no ~ q<fr/T e ~ y/j/u^/k±L. 



n 



(31) 
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FIG. 10. (Color Online) Eigenfunctions of the perturbed potential in the biased plasma config- 
uration (case shown in Fig. [9l bottom) for axial mode number n z = 0,1, azimuthal harmonics 
m e = 3,12,20. 

where 7 and cu* are the growth rate and frequency of the fastest linear mode in the system. 

In the unbiased case (Fig. [91 top), the maximum growth rate is achieved for modes with 
nig ~ 20 — 30, so k± ~ mg/r ~ 100 m _1 . The frequency is close to the growth rate for 
these modes, so ~ 0.5 — 1. The background density gradient scale length near the 

radial location of interest (cathode edge at ~ 28 cm where the equilibrium density gradient 
is mostly localized) is L n = no/n' Q ~ 0.1 m. Therefore, both the simple mixing length 
argument and the weak turbulence theory give a similar estimate for the saturated level 
of turbulence, n/riQ ~ <l4>/T e ~ 10%. This estimate is close to the observed amplitude of 
fluctuations in LAPD measurements and in the nonlinear simulations of LAPD discharg 
The diffusion coefficient estimate based on the mixing length argument, D ~ j/kj_ ~ 2 m 2 /s, 
is close to the value calculated from a saturated state in a self-consistent nonlinear simulation 
D ~ 3 m 2 /s^. This value is comparable to Bohm diffusion, Db ~ 8 m 2 /s, and diffusive 
transport with a Bohm diffusion coefficient has been found to describe the measured profiles 
well in the unbiased case^. 

The biased configuration has reduced radial transport due to strong azimuthal flows and 
is better described by classical diffusion coefficient^. Detailed analysis of this case requires 

22 



a self-consistent nonlinear simulation that takes into account the average radial electric field 
profile; this will be the subject of future work. 

VI. EFFECT OF ION-NEUTRAL COLLISIONS 

The results presented in previous sections do not include ion-neutral collision terms that 
enter the vorticity equation Eq. OH]). The general effect of the v in term is to damp the 
vorticity perturbations (as can be seen from Eq. and to stabilize the wave. In Fig. [TTj 
variation of the real frequency and growth rates for the three modes (drift, KH, and IC) is 
shown as a function of the ion-neutral collisionality parameter u in . Each branch is taken at 
a fixed azimuthal mode number mg that corresponds to the maximum growth rate without 
neutrals (same solution as in Fig. [S]), except for the interchange branch, where mg = 10 is 
chosen) . All of the frequencies and growth rates are normalized to the corresponding values 
at v in = 0. 

When ion-neutral collisions are included, the drift wave growth rate decreases and the 
mode can be completely stabilized at sufficiently high neutral density. For a typical LAPD 
discharge, the rough estimate of is ~ 5 x 10 11 cm -3 — which translates into vi n ~ 
2 x 10~ 3 fi.;. At these values of n^, the effect of the neutrals on the linear stability is 
relatively weak. To completely stabilize the drift mode, should be larger by a factor of 
10 (Fig. HH red). However, due to significant uncertainty in the values of neutral density 
in LAPD, ion- neutral collisions can potentially be important. More importantly, initial 
nonlinear simulations using BOUT show that even at the values near the estimated z/j n ~ 
2 x 10 -3 fi;, the neutral damping is important for the dynamics of the self-generated zonal 
flows^. 

Compared to the drift mode, the KH instability is more strongly affected by the ion- 
neutral collisions. Compared to the neutral-free case, at the estimated for LAPD level of 
ion-neutral collisions, the growth rate drops by ~35% and the mode is completely stabilized 
at Ufo/Sli ~ 0.006. 

The interchange mode turns out to be weakly affected by ion-neutral collisions. For all 
three instability branches, the frequency of the mode remains nearly constant in the range 
of relevant values of neutral collisionality (Fig. [TTJ, dashed lines). 
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FIG. 11. (Color Online) Effect of ion-neutral collisions on the Kelvin-Helmholtz (KH), drift wave 
(DW) and interchange (IC) branches of the dispersion relation. All scans are normalized to the 
corresponding value at Vi n = 0. Solid line - growth rate, dashed - frequency. Estimate for typical 
LAPD parameters: Vml^i ~ 2 x 1CT 3 . 

VII. CONCLUSIONS 

The 3D initial value fluid code for tokamak edge plasma has been adapted to LAPD 
geometry. A separate eigenvalue solver for BOUT set of linearized equations has been 
developed for an independent verification of BOUT results when an analytic solution is not 
available. Background flow terms have been added to BOUT equations to allow simulation 
of flow-driven instabilities. Periodic boundary conditions has been adopted in the parallel 
direction as a first step. A more realistic model of sheath boundary conditions will be 
implemented in future simulations to capture the effect of the parallel boundary on the 
dynamics of the average radial electric field. 

Starting from a system of 3D plasma fluid equations, the derivation of a dispersion relation 
is presented that includes three plasma instability branches: resistive drift mode, Kelvin- 
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Helmholtz mode, and interchange mode; the latter two driven by plasma azimuthal flow. It 
is demonstrated that for LAPD parameters the growth rates for all three branches may be 
comparable, so all three physical mechanisms are potentially important. Interaction with 
neutrals, for the estimated LAPD neutral density, does not significantly affect the linear 
stability of considered modes. However, neutral dynamics can be important for the zonal 
flow generation in nonlinear simulations. 

The initial value solution obtained with BOUT accurately reproduces analytic calcu- 
lations of the properties of the three instabilities, including growth rates, frequencies and 
eigenfunctions. The code solution is in full agreement with analytic and eigenvalue solutions, 
for both model profiles and experimentally relevant profiles, which lends confidence for pro- 
ceeding with nonlinear simulations and validation of BOUT against LAPD measurements. 
Aspects of these linear theoretical estimates (dominant mode numbers, mode frequency, 
and quasilinear estimates of fluctuation amplitudes and diffusion coefficient) are consistent 
with the experimental measurements in LAPD. However, more detailed comparison with 
experiment requires self-consistent nonlinear simulations and this work is underway. Initial 
nonlinear calculations based on the model discussed in this work, and detailed comparisons 
with experimental data, will be presented in a companion paper—. 
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Appendix A: Derivation of fluid equations 

The perpendicular component of the current in Eq. (jlj) is found from the fluid equation 
for the ions Eq. ([3]). Note that the viscosity tensor n and ion pressure terms are dropped 
here, since we neglect the ion temperature effects in this work. Solving it for ion velocity Vj 
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in the Lorentz term, we obtain Vjj_ as a sum of the E x B, polarization and Pedersen drifts: 

V;_L = V^ + Vpj + V/j (Al) 

where v E = cE x B/B 2 , v pi = — (d t + v< • V) v<, v fi = j^-B x Vi/B. 

The main contributions to the perpendicular part of the current divergence come from 
the ion polarization current (the electron olarization drift is smaller by mass ratio) and 
Pedersen current: 

V ■ Ji V ■ (enVpj + envfi) 

= ■ < enb x I — + v* ■ VVj j + enu in b x v, ( 

~ ^- V • {enb x (d t + w E ■ V) v E + enz/ in b x v £ } 

To make the linear expansion of the current continuity equation exactly equivalent to 
the linearized BOUT vorticity equation discussed below, we employ the approximation V ■ 
(nVj) ~ \e ■ Vn (well satisfied for typical LAPD parameters), the same way it is done in 
previous work by Simakov and Cattc3 (Eq. D3): 

9 

V ■ J ± « -- ■ {{d t + v £ • V) (nV ± 0) + i/ in nV_L0} (A2) 

Substituting this expression in the charge conservation equation Eq. combining with 
the continuity equation Eq. ([T]) and parallel projection of the electron momentum equation 
Eq. and linearizing, we obtain: 

d t N + b x V±0 O • ViV = -b x V±0 • ViVo - iV V||i;||e 

T 

d t V\\ e + b X V ± O ■ VV|| e = -/i-r^V||iV + /Lt V|| - V e V\\ e (A3) 

iV 

^oV[|V||e = -V± • (N d t V x <P + d t NV ± fo 

+b x V ± O ■ V (iV o V ± o ) 
+b x V ± 0o ■ V (iV o V ± 0) 
+b x V ± • V (iV o V ± o ) 
+b x V ± O ■ V (iVV ± o ) 

+AT O f in V±0o + N v in V±(j) + Nv in ¥ ± 4 

We project these equations on cylindrical coordinates (r, 9, z) and assume the fluctuations 
are of the form /(x) = /(r) exp(img9 + ik\\z — icut). Solving the first two equations for N and 
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U|| e , and substituting them in the current equation, we obtain ID equation for the perturbed 
potential: 

C 2 (r)0" + C^r)<f>' + C o (r)0 = 0, (A4) 



C 2 (r) = {v in - iu) (A5) 



Ci(r) = {v in - iu) ( + ) + im e -^—4>' (A6) 



C (r) = iv in -iu] 



where 



("^ + Wo (; - ^) + (A^'o)') (A7) 
- r^g - r 2 C - ^ Wo)' + (A8) 
-Hfc[| A„ + ime-Ajv0o0o, (A9) 



1 _ ZeO ™fl 

A„(r,u)) = ik\\n .^ L " 2 Te0 



Ajv(r,o;) - 



u [y e — iu) + ikjnT e0 

T N ~ m «l/ 

£n = — 7w, W = W 0o 

Aq r 



(A10) 



Appendix B: Derivation of the vorticity equation 

Expanding the charge conservation equation V ■ J = as described in section [Til we can 
write 

= V • J 1 1 + V • J± = V||(JVu||) - V± • {(d t + w E ■ V) (nV ± 0) + v in nV±</>} (Bl) 

Introducing the potential vorticity defined as w = Vj_ ■ (AVj_0), we can rewrite the second 
term: 

-V± • {(ft + v E • V) (AV ± 0) + ^„nV ± 0} 

= — d t zu — \ E ■ Vw — Vj_Vb : Vj_(AVj_0) — v in w 

= -d t zu - \ E ■ Vw - Vxy E ■ V±NV±(p - NV±v E ■ V ± V ± - v in w (B2) 
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The fourth term in this expression vanishes: 



V±v B : V±V±0 = Vj_(b x V±0) : V±V±0 
= I (Vi(V ± ■ b x V ±< f>) - (VlV ± 0) ■ (b x V ± 0) - Vi(b x V ± 0) ■ V ± 0) 

= ~ (-(VlV ± 0) • (b x V ± 0) - (b x V ± Vi0) • V ± 0) = o 



(B3) 



The third term in Eq. (1B2I) can be simplified as follows: 



V±V£ : V±iVV ± = {V±0 • V±(b x V±0)} • V±N 



(V±N x b) • (V±0 • V ± V ± 0) 

^(ViV x b) • (V ± V ± 2 ) = ^(ViV x b) • V ± v| 



(B4) 



Collecting all terms, we can write the equation for the evolution of potential vorticity: 



Appendix C: Parameters and profiles for the benchmark case 

Parameters and profiles used for the simulation are presented in Fig. [HJ 

Common parameters for all 3 cases (drift wave, Kelvin-Helmholtz, interchange): 

Helium plasma, once ionized Z = 1 

Radial interval r a < r < r&, r a = 0.15 m, r& = 0.45 m 



{d} = {-5.4638, 124.624, -882.24, 2863.636, -4436.36, 2666.664}, n = 2.5 x 10 18 m" 3 . 

Different parameters for each of the 3 cases: 
Drift wave case: n z = 0.5, 0o( r ) — 0. 
Kelvin-Helmholtz case: n z = 0, 



d t w = - 




B Q = 0.04 T, T e = 5 eV, v in = 0, L z = 17 m 

Density profile is a polynomial fit to the experimental profile rii(r) = no Yll 



Mr) = <Px {F{x - 1) + F(x + 1) - 2F(x)) , 



(CI) 
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x = 4(r - r a )/(r b - r a ) - 2, w = 0.8, <\> x = 50 V. 

fr\ 2 

Interchange case: n z = 0, 4>o( r ) = 4>x [ — , 4>x = 50 V. 

Boundary conditions: periodic in the azimuthal and axial directions; <f>(r a ) = 0(r&) = 
radially. 
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